#include "gbn_cuda.cuh"
using namespace std;


// ------------------------ CUDA设备端数学函数宏定义 --------------------------
/**
 * CUDA指数函数适配：解决float/double类型函数名差异问题
 * @param x：输入指数（支持Float模板类型，如float/double）
 * @return 指数计算结果（float类型调用expf，double类型调用exp，保证精度与性能）
 * @note 1. 依赖C++17 is_same_v编译期类型判断；2. 避免手动分支判断，减少代码冗余
 */
#define CUDA_EXP(x) (is_same_v<Float, float> ? expf(x) : exp(x))



// ------------------------ 依赖的CUDA设备端辅助函数声明 --------------------------
template <typename Float>
__device__ Float gbn0_gaussian(Float distance_sq, Float ss2InvNeg);
template <typename Float>
__device__ Float gbn0_toroidalDistance(Float delta);
template <typename Float>
__global__ void initializeGradientsKernel(Float* d_gradients, int num_points, int dimensions);
template <typename Float>
__global__ void computeGradientKernel(const Float* d_points, Float* d_gradients, 
                                     int num_points, int dimensions, 
                                     Float ss2InvNeg, bool is_toroidal);
template <typename Float>
__global__ void updatePointsKernel(Float* d_points, const Float* d_gradients, 
                                   int num_points, int dimensions, 
                                   Float scale, bool is_toroidal);
template <typename Float>
__global__ void boundedDomainForceKernel(const Float* d_points, Float* d_gradients, 
                                         int num_points, int dimensions, Float sigma);

// ------------------------ 依赖的CUDA设备端辅助函数实现 --------------------------
 /**
 * 高斯函数值计算：基于预计算参数ss2InvNeg，减少设备端重复计算
 * @param distance_sq：点对的平方距离（各维度距离平方和，避免sqrt操作，提升性能）
 * @param ss2InvNeg：高斯核预计算参数（= -1/(2*sigma²)，由主机端传入，减少设备端计算量）
 * @return 高斯权重值（范围0 < 权重 ≤ 1，近邻点权重高，远邻点权重低）
 * @note 1. 支持任意维度场景（输入为累积平方距离）；2. 完全对齐gbn0的高斯核计算逻辑
 */
template <typename Float>
__device__ Float gbn0_gaussian(Float distance_sq, Float ss2InvNeg) {
    return CUDA_EXP(ss2InvNeg * distance_sq); // 直接使用预计算参数，减少设备端计算量
}

// 设备端辅助函数：处理环形域距离(适用于单周期范围)
/**
 * 环形域距离修正：对齐gbn0的边界逻辑，支持任意维度
 * @param delta：某维度的坐标差（点i - 点j）
 * @return 修正后的最短环形距离（含符号，匹配gbn0方向特性）
 */
template <typename Float>
__device__ Float gbn0_toroidalDistance(Float delta) {
    // 完全复刻gbn0的环形域修正逻辑：处理±0.5边界
    if (delta < -static_cast<Float>(GBN0_HALF)) delta += static_cast<Float>(GBN0_ONE);
    if (delta >= static_cast<Float>(GBN0_HALF)) delta -= static_cast<Float>(GBN0_ONE);
    return delta;
}

// 设备端核函数：初始化GPU端梯度数组为零
/**
 * 梯度数组清零：为每次迭代的梯度计算提供初始状态
 * @param d_gradients：输出，GPU端梯度数组指针（长度=num_points×dimensions，存储每个点的各维度梯度）
 * @param num_points：输入，点集总数量（如4096、8192）
 * @param dimensions：输入，点集维度（如2D、3D，最大支持MAX_DIMENSIONS）
 * @note 1. 线程索引覆盖所有梯度元素（避免未初始化的脏数据）；2. 与computeGradientKernel配套使用
 */
template <typename Float>
__global__ void initializeGradientsKernel(Float* d_gradients, int num_points, int dimensions) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < num_points * dimensions) {
        d_gradients[idx] = static_cast<Float>(0.0);
    }
}

// 设备端核函数：计算点集的梯度（GBN算法核心逻辑）
/**
 * 梯度计算：基于高斯核加权，实现点集的分散化优化
 * @param d_points：输入，GPU端当前点集数组（长度=num_points×dimensions，存储每个点的坐标）
 * @param d_gradients：输出，GPU端梯度数组（长度与d_points一致，存储计算后的各点梯度）
 * @param num_points：输入，点集总数量（N）
 * @param dimensions：输入，点集维度（D）
 * @param ss2InvNeg：输入，高斯核预计算参数（= -1/(2*sigma²)，对齐gbn0逻辑）
 * @param is_toroidal：输入，边界模式开关（true=环形域，false=有界域）
 * @note 1. 单线程处理单个点的梯度计算（避免线程竞争）；2. 环形域场景调用gbn0_toroidalDistance修正距离；3. 自动跳过点自身（i==j）的无效计算
 */
template <typename Float>
__global__ void computeGradientKernel(
    const Float* d_points,        // 输入：GPU端点集数组（长度=num_points*dimensions）
    Float* d_gradients,           // 输出：GPU端梯度数组（与d_points同长度，存储每个点的梯度）
    int num_points,               // 输入：点集总数量（N）
    int dimensions,               // 输入：点集维度（如2D、3D，最大支持32维）
    Float ss2InvNeg,              // 输入：高斯预计算参数（= -1/(2*sigma²)，对齐gbn0）
    bool is_toroidal              // 输入：是否启用环形域（true=环形，false=有界域）
) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= num_points) return;

    // 初始化当前点的梯度
    Float grad[32];  // 栈上数组，存储点i的各维度的梯度（最大支持32维，适配高维场景）
    for (int d = 0; d < dimensions; d++) {
        grad[d] = static_cast<Float>(0.0);  // 所有维度梯度初始化为0
    }


    // 获取当前点坐标
    const Float* x_i = &d_points[i * dimensions];   //找到第i个点在d_points数组中的第一个坐标元素的地址

    // 遍历所有其他点
    for (int j = 0; j < num_points; j++) {
        if (i == j) continue;   // 跳过自身点，避免自贡献

        // 获取对比点坐标
        const Float* x_j = &d_points[j * dimensions];   //找到第j个点在d_points数组中的第一个坐标元素的地址

        // 计算各维度距离
        Float dist_sq = static_cast<Float>(0.0);   // 点i与j的平方距离（避免sqrt，减少计算量）
        Float delta[32];        // 存储各维度修正后的坐标差
        for (int d = 0; d < dimensions; d++) {
            delta[d] = x_i[d] - x_j[d]; // 原始坐标差（点i - 点j）
            // 环形域修正：调用gbn0风格的距离函数；有界域不修正（保留原始差值）
            if (is_toroidal) {
                delta[d] = gbn0_toroidalDistance(delta[d]);
            }
            dist_sq += delta[d] * delta[d]; // 累积平方距离
        }

        // 计算高斯权重（近邻点权重高，远邻点权重低）
        Float g = gbn0_gaussian(dist_sq, ss2InvNeg);  // 调用设备端辅助函数，得到权重值（0 < g ≤ 1）

        // 累积梯度：加权坐标差贡献到各维度梯度
        for (int d = 0; d < dimensions; d++) {
            grad[d] += g * delta[d];  // 权重×坐标差 = 该对比点j对i的梯度贡献
        }
    }

    // 应用梯度公式系数
    //Float coeff = -PI / (Float)num_points;  // 梯度系数（来自GBN算法理论，用于归一化梯度大小）
    for (int d = 0; d < dimensions; d++) {
        // 计算最终梯度值，写入全局数组（供后续坐标更新使用）
        d_gradients[i * dimensions + d] = grad[d];
    }
}

// 设备端核函数：根据梯度更新GPU端点集坐标
/**
 * 坐标更新：基于梯度下降法，调整点集位置以优化分布
 * @param d_points：输入输出，GPU端当前点集数组（需修改，存储更新后的坐标）
 * @param d_gradients：输入，GPU端梯度数组（computeGradientKernel的输出结果，只读）
 * @param num_points：输入，点集总数量
 * @param dimensions：输入，点集维度
 * @param scale：输入，学习率（控制单次更新步长，防止优化震荡或收敛过慢）
 * @param is_toroidal：输入，边界模式开关（true=环形域，false=有界域）
 * @note 1. 环形域场景对齐gbn0的[0,1)边界修正逻辑；2. 有界域场景采用[0,1]范围裁剪；3. 单线程处理单个点的所有维度更新
 */
template <typename Float>
__global__ void updatePointsKernel(
    Float* d_points,        // 输入输出：GPU端当前点集数组（需修改，无const）
    const Float* d_gradients, // 输入：GPU端梯度数组（只读，不修改）
    int num_points,         // 输入：点集总数量（如4096）
    int dimensions,         // 输入：点集维度（如2D、3D）
    Float scale,               // 输入：学习率（控制单次更新的步长）
    bool is_toroidal        // 输入：边界模式（true=环形域，false=有界域）
){
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= num_points) return;

    for (int d = 0; d < dimensions; d++) {// 遍历当前点的所有维度（如2D的x、y）
        int idx = i * dimensions + d;       // 计算当前维度坐标在数组中的索引（点i的第d维）
        // 核心更新公式：新坐标 = 旧坐标 + 学习率 × 梯度值
        Float new_coord = d_points[idx] + scale * d_gradients[idx];

        // 边界处理：保留环形域/有界域双接口
        if (is_toroidal) {
            // 环形域：复刻gbn0的边界修正逻辑（处理[0,1)范围）
            if (new_coord < GBN0_ZERO) new_coord += GBN0_ONE;
            if (new_coord >= GBN0_ONE) new_coord -= GBN0_ONE;
        } else {
            // 有界域：保留原代码的[0,1]裁剪逻辑，不修改
            new_coord = fmaxf(GBN0_ZERO, fminf(GBN0_ONE, new_coord));
        }
        // 写入更新后的坐标
        d_points[idx] = new_coord;
    }
}

// 设备端核函数：为有界域场景添加边界引力（防止点聚集在边界）
/**
 * 边界引力计算：补充有界域的边界约束，优化点集边缘分布
 * @param d_points：输入，GPU端当前点集数组（存储各点坐标，只读）
 * @param d_gradients：输入输出，GPU端梯度数组（在原有梯度基础上累加边界引力）
 * @param num_points：输入，点集总数量
 * @param dimensions：输入，点集维度
 * @param sigma：输入，高斯核带宽（与computeGradientKernel的sigma一致，保证引力与梯度的尺度匹配）
 * @note 1. 仅在is_toroidal=false（有界域）时调用；2. 基于erf函数计算边界引力，靠近边界时产生反向拉力；3. 引力强度通过固定系数（0.1）调节，避免影响主梯度
 */
template <typename Float>
__global__ void boundedDomainForceKernel(const Float* d_points, Float* d_gradients, 
                                         int num_points, int dimensions, Float sigma) 
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= num_points) return;

    const Float* x_i = &d_points[i * dimensions];
    Float sigma_half = sigma * static_cast<Float>(0.5);

    for (int d = 0; d < dimensions; d++) {
        Float x = x_i[d];
        // 计算边界引力（基于erf函数）
        Float term1 = CUDA_EXP(-powf((static_cast<Float>(1.0) - x) / (static_cast<Float>(2.0) * sigma), 2.0));
        Float term2 = CUDA_EXP(-powf(x / (static_cast<Float>(2.0) * sigma), 2.0));
        Float force = term1 - term2;

        // 累加到梯度
        d_gradients[i * dimensions + d] += force * static_cast<Float>(0.1);  // 调整引力强度
    }
}


// -------------------------- GBNPointSet 类成员函数实现 --------------------------
// 类构造函数：初始化点集规模、域范围与GPU资源
/**
 * 构造函数：完成基类初始化、维度合法性校验与GPU梯度内存分配
 * @param num_points：输入，点集总数量（如4096、8192）
 * @param domain_ranges：输入，各维度的域范围（如2D场景{1.0f, 1.0f}表示x/y轴均为[0,1)）
 * @throw invalid_argument：维度超出1~MAX_DIMENSIONS范围时抛出
 * @note 1. 调用基类GBNBasePointSet的构造函数初始化点集核心参数；2. 自动调用initDeviceGradients分配GPU梯度内存；3. 初始化快照默认参数（如输出路径、格式）
 */
template <typename Float>
GBNPointSet<Float>::GBNPointSet(int num_points, const vector<Float>& domain_ranges)
    : GBNBasePointSet<Float>(num_points, domain_ranges),  // 调用基类构造（初始化点集规模与域范围）
      d_gradients_(nullptr) {
    // 校验维度合法性（与MAX_DIMENSIONS一致，避免GPU栈溢出）
    if (this->dimensions <= 0 || this->dimensions > MAX_DIMENSIONS) {
        throw invalid_argument(
            "Dimensions must be between 1 and " + to_string(MAX_DIMENSIONS) + 
            " (current: " + to_string(this->dimensions) + ")"
        );
    }

    // 初始化GPU梯度内存与快照默认参数
    initDeviceGradients();
}

// 类析构函数：释放GPU梯度资源
/**
 * 析构函数：安全释放类自身管理的GPU资源，避免内存泄漏
 * @note 1. 仅释放d_gradients_（GPU梯度数组），基类析构函数自动释放d_points、d_points_ref等GPU资源；2. 调用freeDeviceGradients确保释放逻辑安全（空指针检查、释放后置空）
 */
template <typename Float>
GBNPointSet<Float>::~GBNPointSet() {
    // 释放GPU梯度资源（基类析构会自动释放d_points/d_points_ref等）
    freeDeviceGradients();
}

// 类成员函数：GBN算法核心优化接口（无参数，复用类成员变量）
/**
 * 优化主函数：实现完整的梯度下降迭代流程，输出优化后的点集
 * @throw invalid_argument：优化参数（sigma、学习率、迭代次数）非法时抛出
 * @throw runtime_error：CUDA核函数执行失败或设备同步错误时抛出
 * @note 1. 前置调用validateOptimizeInput校验参数合法性；2. 自动完成CPU→GPU的初始点集复制；
 *       3. 迭代过程包含梯度初始化、计算、更新、学习率衰减；4. 支持快照生成（按视觉变化触发）；
 *       5. 迭代结束后自动将GPU端优化结果复制回CPU（h_points）
 */
template <typename Float>
void GBNPointSet<Float>::optimize() {
    // 步骤1：校验优化参数合法性（sigma、学习率、迭代次数）
    validateOptimizeInput();

    // 步骤2：将主机端初始点集复制到GPU（复用基类数据传输接口）
    this->copyHostToDevice();          // CPU→GPU（d_points）
    this->copyDeviceToRef();           // GPU→参考点集（d_points_ref，快照用）
    int total_size = this->num_points * this->dimensions;

    // 快照前置校验：仅2D支持可视化（TXT支持任意维度）
    int current_format = SNAPSHOT_FORMAT_MASK;
    if ((current_format & SNAP_FORMAT_PNG) && this->dimensions != 2) {
        cout << "[Warning] PNG visualization only supports 2D! Will disable PNG output." << endl;
        // 修正格式掩码：移除PNG，仅保留TXT
        current_format &= ~SNAP_FORMAT_PNG;
        // 3. 将修改后的格式同步到基类（关键！确保快照按新格式输出）
        this->setFrameFormat(current_format);
    }

    // 设备同步与错误检查（复用基类接口）（确保初始化阶段无错误）
    cudaDeviceSynchronize();
    this->checkCudaErr("CUDA error during initialization (before iterations)"); 

    // 步骤3：设置CUDA核启动参数（块大小固定为CUDA_BLOCK_SIZE，网格大小自适应）
    dim3 block(CUDA_BLOCK_SIZE);
    dim3 grid_points((this->num_points + block.x - 1) / block.x);       // 点级核（如computeGradient）的网格
    dim3 grid_total((total_size + block.x - 1) / block.x);              // 梯度级核（如initializeGradients）的网格


    // 高斯预计算参数：ss2InvNeg = -1/(2*sigma²)
    Float ss2InvNeg = -1.0f * this->num_points / (2.0f * this->sigma * this->sigma);
    Float scale = this->learning_rate;

    // 生成快照（复用基类的 generateSnapshot 接口）
    this->generateSnapshot(0);


    // 梯度下降迭代
    for (int iter = 0; iter < this->max_iterations; iter++) {
        // 1. 初始化梯度
        initializeGradientsKernel<<<grid_total, block>>>(d_gradients_, this->num_points, this->dimensions);
        this->checkCudaErr("initializeGradientsKernel failed at iteration " + to_string(iter));

        // 2. 计算点集核函数梯度
        computeGradientKernel<<<grid_points, block>>>(
            this->d_points, d_gradients_,
            this->num_points, this->dimensions,
            ss2InvNeg, this->is_toroidal
        );

        // 3. 有界域边界引力额外处理
        if (!this->is_toroidal) {
            boundedDomainForceKernel<<<grid_points, block>>>(
                this->d_points, d_gradients_,
                this->num_points, this->dimensions,
                sigma
            );
        }

        // 4. 更新点坐标
        updatePointsKernel<<<grid_points, block>>>(
            this->d_points, d_gradients_,
            this->num_points, this->dimensions,
            scale, this->is_toroidal
        );
        this->checkCudaErr("updatePointsKernel failed at iteration " + to_string(iter));

        // 5. 学习率衰减（每1000迭代衰减10%）
        if ((iter + 1) % 1000 == 0) {
            scale *= 0.9f;
            cout << "Iteration: " << iter+1 << ", Current scale: " << scale << endl;
        }

        // 5.6 生成快照（复用基类generateSnapshot接口，按视觉变化触发）
        this->generateSnapshot(iter + 1);

        // 5.7 设备同步与错误检查（复用基类checkCudaErr）
        cudaDeviceSynchronize();
        this->checkCudaErr("CUDA error after iteration " + to_string(iter));
    }

    // 复制结果到主机端
    this->copyDeviceToHost();          // GPU→CPU（h_points）

    // 快照收尾：输出统计信息
    cout << "GBN Optimization Finished!" << endl;
    cout << "Total Iterations: " << this->max_iterations << endl;
    cout << "Total Snapshots Recorded: " << this->recorded_frames << endl;
    cout << "Optimized Points Saved to Host (h_points)!" << endl;
}

// 类成员函数：手动释放GPU梯度内存（安全接口）
/**
 * @brief 手动释放GPU梯度内存（安全释放，避免野指针与重复释放）
 * @note 1. 仅释放类自身管理的 `d_gradients_`（基类GPU资源由基类析构函数释放）
 *       2. 释放后将指针置空，防止后续误操作（如重复释放、访问野指针）
 *       3. 复用基类 `checkCudaErr` 接口，统一错误处理风格
 */
template <typename Float>
void GBNPointSet<Float>::freeDeviceGradients() {
    // 1. 跳过空指针：若梯度数组未分配或已释放，直接返回（避免 `cudaFree` 空指针报错）
    if (d_gradients_ == nullptr) {
        std::cout << "[Info] GPU gradient array (d_gradients_) is already null, skip free." << std::endl;
        return;
    }

    // 2. 执行CUDA内存释放（释放 `d_gradients_` 指向的GPU内存）
    cudaError_t err = cudaFree(d_gradients_);
    // 3. 错误检查：复用基类 `checkCudaErr` 接口，统一抛出异常（确保错误可捕获）
    this->checkCudaErr("freeDeviceGradients: failed to free GPU gradient array");

    // 4. 释放后置空指针：避免后续误操作（如再次调用 `free` 或访问已释放内存）
    d_gradients_ = nullptr;

    // 5. 日志提示：便于调试时确认资源释放状态
    std::cout << "[Info] GPU gradient array (d_gradients_) freed successfully." << std::endl;
}


// 类成员函数：初始化GPU梯度内存（私有辅助函数）
/**
 * 梯度内存分配：为GPU端梯度数组分配内存，适配点集规模
 * @throw runtime_error：cudaMalloc分配GPU内存失败时抛出
 * @note 1. 私有函数，仅在构造函数中自动调用；2. 内存分配大小；
 *       3. 复用基类checkCudaErr检查分配结果
 */
template <typename Float>
void GBNPointSet<Float>::initDeviceGradients() {
    int total_size = this->num_points * this->dimensions;
    cudaMalloc(&d_gradients_, total_size * sizeof(Float));
    this->checkCudaErr("Failed to allocate GPU gradient array (d_gradients)");
}

/**
 * @brief 输入参数校验（优化前前置检查，保障核心逻辑安全）
 * @note 1. 校验范围：基类点集合法性 + 类成员优化参数（sigma/learning_rate/max_iterations）
 *       2. 复用基类成员（h_points/dimensions），无需外部传参，符合类封装原则
 *       3. 非法时抛出带具体原因的 invalid_argument，便于调试定位
 */
template <typename Float>
void GBNPointSet<Float>::validateOptimizeInput() const {
    // -------------------------- 1. 校验基类主机端点集（h_points）合法性 --------------------------
    // 1.1 点集非空检查：未初始化点集（如外部未调用基类 setPoints 接口）直接优化无意义
    if (this->h_points.empty()) {
        throw invalid_argument(
            "Host point set (h_points) is empty! "
            "Please initialize it first via base class interface (e.g., setPoints())."
        );
    }

    // 1.2 点集尺寸与维度匹配检查：h_points的size = 点数量 
    // 例：2D场景100个点 → h_points 长度需为 100×2=200
    //size_t expected_total_size = static_cast<size_t>(this->num_points) * this->dimensions;
    if (this->h_points.size() != static_cast<size_t>(this->num_points)) {
        throw invalid_argument(
            "Host point set size mismatch! "
            "Expected number of points: " + to_string(this->num_points) + 
            ", Actual number of points: " + to_string(this->h_points.size())
        );
    }

    // 1.3 点坐标范围合法性（可选，基于环形域/有界域场景补充）
    // 若基类支持边界模式（is_toroidal），可额外校验坐标是否在 [0, 1) 范围内（匹配 GBN0 标准）
    //遍历每个点的每个维度，通过 Point 的 operator[] 访问坐标
    if (this->is_toroidal) {
        for (size_t i = 0; i < this->h_points.size(); ++i) { // 遍历每个点（Point<Float>）
            const Point<Float>& current_point = this->h_points[i]; // 取当前点的引用（避免拷贝）
            for (int d = 0; d < this->dimensions; ++d) { // 遍历当前点的每个维度
                // 修正：通过 Point::operator[] 访问第 d 维坐标（需确保 Point 类已实现此接口）
                Float coord = current_point[d]; 
                if (coord < GBN0_ZERO || coord >= GBN0_ONE) {
                    throw invalid_argument(
                        "Toroidal domain coordinate out of range! "
                        "Point " + to_string(i) + ", Dimension " + to_string(d) + 
                        ": Expected [0.0, 1.0), Actual: " + to_string(coord) + 
                        "\nHint: Use base class normalizePoints() to adjust coordinates."
                    );
                }
            }
        }
    }

    // -------------------------- 2. 校验类成员优化参数合法性 --------------------------
    // 2.1 高斯核带宽（sigma）：必须为正数（核带宽为负/零 → 高斯函数无衰减，权重计算失效）
    if (sigma <= static_cast<Float>(0.0)) {
        throw invalid_argument(
            "Gaussian kernel bandwidth (sigma) must be positive! "
            "Current value: " + to_string(sigma)
        );
    }

    // 2.2 学习率（learning_rate）：必须为正数（负学习率→梯度反向更新，零学习率→无更新）
    if (learning_rate <= static_cast<Float>(0.0)) {
        throw invalid_argument(
            "Learning rate must be positive! "
            "Current value: " + to_string(learning_rate)
        );
    }
    // 额外约束：学习率不宜过大（避免优化震荡，基于工程经验设上限10.0）
    if (learning_rate > static_cast<Float>(10.0)) {
        throw invalid_argument(
            "Learning rate is too large (may cause optimization divergence)! "
            "Current value: " + to_string(learning_rate) + ", Recommended max: 10.0"
        );
    }

    // 2.3 最大迭代次数（max_iterations）：必须为正整数（零/负迭代→优化逻辑无法执行）
    if (max_iterations <= 0) {
        throw invalid_argument(
            "Max iterations must be a positive integer! "
            "Current value: " + to_string(max_iterations)
        );
    }
    // 额外约束：迭代次数不宜超过硬件承载上限（避免长时间占用GPU，参考 MAX_ITERATIONS 常量）
    if (max_iterations > MAX_ITERATIONS) {
        throw invalid_argument(
            "Max iterations exceeds hardware limit! "
            "Current value: " + to_string(max_iterations) + ", Limit: " + to_string(MAX_ITERATIONS) + 
            "\nHint: Modify MAX_ITERATIONS constant if you need more iterations."
        );
    }


    // -------------------------- 3. 校验基类GPU资源初始化状态（可选，防重复分配） --------------------------
    // 若优化函数中需复用基类GPU内存（如 d_points），可提前校验是否已分配（避免重复调用 copyHostToDevice）
    if (this->d_points == nullptr) {
        throw invalid_argument(
            "Base class GPU point array (d_points) is uninitialized! "
            "Please call copyHostToDevice() before optimize()."
        );
    }
}

// 关键：显式实例化 float 版本的 GBNPointSet
template class GBNPointSet<float>; // 生成 GBNPointSet<float> 的所有成员函数代码（包括 optimize()）

// 若后续需要 double 版本，可添加：
// template class GBNPointSet<double>;